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solutions of the time-dependent Schrodinger equation. The generalization yields numerical solutions 
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employ r = M = 1. We note dramatic improvement in the attainable precision (circa 10 or greater 
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short-range interactions, wavepacket scattering, and long-time studies of decaying systems. 
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I. INTRODUCTION 

Whereas there are a number of examples of exact an- 
alytic solutions of time-independent problems in quan- 
tum mechanics, such solutions of time-dependent prob- 
lems are few. The analytic solutions of both types that 
do exist tend to provide approximate models to actual 
physical systems. No doubt such solutions are instruc- 
tive for gaining insight into the behavior of the physi- 
cal systems that they describe. Nevertheless because the 
models themselves are often approximations and because 
one wishes to describe real systems as precisely as possi- 
ble, one also relies on accurate numerical methods. 

The time-dependence of nonrelativistic quantum sys- 
tems, which is the focus of this paper, has become 
important in diverse areas of atomic and subatomic 
physics. Examples of these include the study of nu- 
clear processes such as the decay of unstable nuclei and 
associated phenomena like atomic ionization [J, [2J and 
bremsstrahlung [!, 0, [f| , the study of fundamental pro- 
cesses necessary for quantum computing [f| , the study of 
mesoscopic physics or nanophysics devices Q, and the 
motion of atoms in a trap. A reliable and accurate nu- 
merical determination of the time-dependent wave func- 
tion such as we discuss in this paper will no doubt be 
necessary and/or helpful in making advances in the un- 
derstanding of a variety of quantum processes. 

In this paper we consider the numerical solution of 
the time-dependent nonrelativistic Schrodinger equation. 
Much has been learned about basic scattering processes 
from the numerically generated solutions of traveling 
wavepackets as they pass through a potential region Q, 
as well as the time-evolution of unstable quantum pro- 
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cesses However, the methods used in the past, and 
still employed currently, are limited in that the solutions 
often degrade after a certain time interval, so that they 
reduce to noise. Furthermore for processes in which the 
wave function spreads or travels away from the source 
one often requires such a large number of space steps 
that the computation becomes prohibitive. 

The goal of this study is to improve the existing stan- 
dard approach by allowing for relatively large step sizes 
both in time and in space and thus to reduce the number 
of basic arithmetical calculations while obtaining more 
accurate solutions. We have been able to make signifi- 
cant improvement to one conventional approach, viz., the 
Crank-Nicolson (CN) implicit integration scheme for the 
time-dependent Schrodinger equation. Many years ago 
the CN approach was shown to be successful in the study 
of wavepacket scattering in one dimension by Goldberg 
et al. :8j. In recent years the CN method continues to 
be employed for its space- and/or time-development al- 
gorithm to study various time-dependent problems. See, 
for example, Refs. [l(| EH E2, EH- The attractive as- 
pect of this method is that the solution is constrained to 
be unitary at every time step. It is this constraint that 
makes the solution stable regardless of time- or space-step 
size. Although the evolution of the solution is unitary, 
the wave function is not correct if the step sizes are too 
large. The error is of 0((Axf , (At) 3 ), where Ax and At 
are the spatial and temporal step size, respectively. 

The method was successfully generalized to two di- 
mensional scattering by Galbraith et al. 1J| and, more 
recently, to multichannel scattering (l5j . Furthermore al- 
ternative methods which are fast computationally were 
introduced by Kosloff and Kosloff [16]. These involve the 
fast Fourier transform of the kinetic energy operator of 
the Schrodinger e quat ion. Variants of this method were 
discussed in Ref. Although this approach is fast 

and is able to handle large time intervals in one pass, it 
is not unitary and does require a large number of space 
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intervals. 

Improved CN algorithms have been discussed by a few 
authors. Mi§ucu et al. [f| introduced a seven-point for- 
mula for the second-order spatial derivative with error 
of 0((Ax) 6 ) and an improved time-integration scheme 
with an error of 0((At) 5 ). They claim to obtain two 
orders of magnitude improvement in the time advance. 
Moyer [l8[ uses a Numerov scheme for the spatial inte- 
gration method with error of 0((Ax) 6 ) but has a one- 
stage time evolution giving the CN precision in time, 
i.e., with an error of 0((At) 3 ). Moyer also introduces 
transparent boundary conditions for unconfmed systems. 
We have found these very useful for making long-time or 
large-space problems tractable Q , but we will not discuss 
such boundary conditions further in this paper. Puzynin 
et al. [l9|, H(| indicate how to generalize the time devel- 
opment to higher order, but do not discuss spatial inte- 
gration. 

In this paper we present a generalization of the often- 
used CN method of obtaining numerical solutions of 
the time-dependent Schrodingcr equation. The gener- 
alization yields numerical solutions accurate to order 
(Ax) 2r_1 in space and (At) 2M in time for any positive 
integers r and M, while CN employ r = M = 1. By 
appropriate choice of r and M the improvement can be 
of such a nature that hitherto computationally unfeasible 
problems become doable, and solutions with low to mod- 
est precision can now be obtained extremely accurately. 

In the following we consider the generalization of spa- 
tial integration in Sec. [II] and the generalization of the 
time integration in Sec. IIIII In Sec. IIVI we discuss errors 
and a way of dealing with a particular type of bound- 
ary condition. We study specific examples to illustrate 
the improvement of the generalizations over the standard 
CN procedure in Sec. [V] Some general observations and 
conclusions are made in Sec. ED 
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and (j>(x) is a given wave function at initial time t . In 
this section we use the standard time-advance procedure 
of the CN method, but generalize the spatial integration. 
In Sec. |TIT] we generalize the time-evolution procedure. 

The time evolution of the system can be expressed in 
terms of an operator acting on the wave function at time 
t which gives the wave function at a later time t + At 
according to the equation 



ip(x,t + At) 



= -iHAt/h 



il>{x,t). (2.3) 



The time-evolution operator e iHAt/h can ^ e ex p anc Jed 
to give a unitary approximation of the operator by setting 



iHAt/h = 1 - \iHAt/h 
1 + UHAt/h 



+ 0{{Atf). (2.4) 



Inserting the approximate form of the operator into 
Eq. (|2.3|) . we obtain the equation 



+ hiHAt/nj 4>{x,t+At) = (l - ^iHAt/hj tp(x,t), 

(2.5) 

with an error of 0((At) 3 ). Here we focus on the second- 
order spatial derivative in H of Eq. (|2.2[) and leave im- 
provements with respect to the time derivative to Sec. lIIII 
We generalize the usual three-point formula and the 
seven-point formula of Mi§icu et al. 5j, for the second- 
order derivative to a (2r + l)-point formula. Such a for- 
mula has the form 



II. SPATIAL INTEGRATION 

We describe a general procedure for solving the one- 
dimensional time-dependent Schrodinger equation 

i^-flj4t) = 0, VOMo) =<Kx), (2.1) 
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k=r 
h 2 ^ 



c^ ] y{x + kh) + 0{h 2r ), (2.6) 



(r) 

where c k are real constants. To obtain the coefficients 
c£ we make expansions 



y(x + kh) = y(x) + (kh)yW(x) + l(kh) 2 y^(x) + ■■■+ ± (kh) 2r+1 y^ (x) + 0(h 2r + 2 ) 
y(x - kh) = y(x) - (kh)y^(x) + ±(khfyW(x) ■ ■ ■+ { -0^{kh) 2r+l y^ +l \x) + 0(h 2r + 2 ), 

for k = 1, 2, . . . , r; y^ denotes the ith derivative with respect to x. When we add the two equations, the terms with 
odd-order derivatives cancel, resulting in the equation 

2^y (2) (-) + 2^y<%) + • ■ • + 2^V 2r) (x) = y(x + kh) + y(x kh) 2y(x) + 0(h 2 ^ 2 ). (2.7) 
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Thus we obtain the system of r equations in r unknowns, i.e., y^ 2k \x) for k = 1, . . . ,r, 

2^ 2 Hx)+2^y^\x) + --- + 2^y^\ X ) = y(x + h) + y(x - h) - 2y(x) 
2^y^ 2 \x) + 2^y^(x) + --- + 2^ y ^(x) = y(x + 2h) + y(x - 2h) - 2y(x) 

2^y (2) (x) + 2 { ^fy^(x) + ■■■ + 2 { j0^yW(x) = y(x + rh) + y(x - rh) - 2y(x). (2. 



We solve these equations to obtain y^ (x). It is evident 
from the terms on the right side of Eqs. (|2.8[) that y^ 2 '(x) 
has the form of Eq. (|2.6[) and the coefficients c k can 



be identified. Because the equations (|2.8[) are invariant 
under the change of h to —h, the coefficients satisfy the 

(r) (r) 

relation c_ k = c k for k = 1,2, ...,r. For example, 
the first seven sets of coefficients (up to the fifteen-point 
formula) are given in Table [T] 
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2 

1925 
7 

3300 



16632 
7 



:»i.s.ss 84084 



TABLE I: The coefficients c[ r) up to r = 7. 

Let us partition the range of x and i values so that 
xj = xq + jAx, j = 0, 1, . . . , J and t n = to + nAt, 
n = 0, 1, . . . , N. The numerical approximation of the 
wave function at a mesh point in space and time is de- 
noted as i\ij yrl « ip(xj,t n ) and we set Vj = V(xj). Using 
expression (12. 6p in Eq. (|2.5[) . we obtain 
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Vjl/>j,n, (2-9) 



for j — to J. The indices in the sums may go out of 
range, so we set ifjj,n = when j < and j > J. Define 



2m(Ax) 2 ' 
and subsequently 



-2 and 



,0) _ 



4)4 r \ (2-10) 
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The notation includes zj 1 ' which is consistent with that 
used in the generalization of the time dependence of the 
wave function discussed in the next section. 

The solution %j)j %n+ i is obtained by solving the system 
of linear equations 

AV n+1 = A*t> n , (2.12) 
where the matrix A is the (2r + l)-diagonal matrix 
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(2.13) 

where the superscript (r> of the au is assumed. The ma- 
trix A* is the complex conjugate of matrix A. The wave 
function at t„+i, i.e., ^f n +x, is a column vector consisting 
of the ipj,n+x as components, and can be determined if 
^ n is known. The matrix equation (|2.12p can be solved 
using standard techniques. 



III. TIME ADVANCE 

In this section we extend the work of Puzynin et 
al. [Hi, Hfj. The basic idea is to replace the exponen- 
tial operator exp(— iH At) by the diagonal Pade approxi- 
mant. The [M/M] Pade approximant of the exponential 
function may be written as 



M 



z _ a + axz + ■ ■ ■ + aMZ 



a " 

m—0 



bo + bxz + 



+ b M z M 



M 



5> 

m'=0 



(3.1) 
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where the a m and the b m > are complex constants. It is 
evident that when z = 0, a /bo = 1, which makes one of 
the coefficients arbitrary. By convention we take bo = 1 
which immediately fixes ao = 1. There arc 2 M constants 
remaining, which can be found from the known coeffi- 
cients of the series expansion of the exponential function, 
giving an error term in Eq. (|3.1[) 0(z 2M+1 ). The prop- 
erty of Pade approximants that can be used to advantage 
is that, if f(z) is unitary, so is its diagonal Pade approx- 
imant [2l|. 

In general we solve for the coefficients a m and b m i by 
multiplying Eq. (|3.ip by the denominator so that 

( E h ™' Zm ' J ( f>^ J = ( E a «^ m J . (3-2) 
\m'=0 / \i=0 / \m=0 / 

where the Cj are known since e z = X)i*lo Z V*" Multi- 
plying out the sums on the left side of Eq. (|3.2p , and 
equating the coefficients of z through z 2M on both sides, 
we obtain 2M equations in 2M unknowns. The last M of 



these equations contain no a m and hence can be solved 
for the b m i, which in turn can be inserted in the first 
M equations to obtain the a rn . The numerator and the 
denominator of the diagonal Pade approximant of the 
exponential function have been studied extensively [2l| . 
When each is factored it is found that the roots of the 
denominator are the negative complex conjugates of the 
roots of the numerator. Thus the [M/M] Pade approxi- 
mant of the exponential function leads to 

M / , (A/)\ 

s=l \1+Z/Z S j 

where z s , s = 1, . . . , M, are the roots of the numera- 
tor, and zi is the complex conjugate of Zs* . These 
roots can be found to a desired precision for virtually any 
value of M. We have found them to 17 digit precision 
for M up to 20, a sample of which for M = 1 to 5, each 
rounded to five decimal places, is given in Table uTl 
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-3.00000 + il.73205 


-3.00000 - il.73205 
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-4.64437 + iO.OOOOO 


-3.67781 - i3.50876 


-3.67781 + i3.50876 






4 


-4.20758 + i5.31484 


-5.79242 + i 1.73447 


-5.79242 - i 1.73446 


-4.20758 - 


i5.31483 


5 


-4.64935 + i7.14205 


-6.70391 +i3.48532 


-7.29348 + iO.00000 


-6.70391 - 


i3.48532 -4.64935 - i7. 14205 



TABLE II: The roots z s of the numerator of the Pade ap- 
proximant of the exponential function for M from 1 to 5. 



We use the Pade approximant to express the time evo- 
lution operator. Define the operator 



1 iHAt/h 

K (M) = Zl 

3 - , iHAt/h 

1 H -(A/) 



(3.4) 



so that 



M 



e -iHAt/h = -Q K (M) + ((At) 2M+1 ). (3.5) 

s=l 

Since * n+ i = e^^At/h^^^ we wr ^ e foe relation 

AI 

*„ +1 =n#i M) *». (3-6) 



8 = 1 



Defining *„+ s / M = Ks * ra +( s -i)/M) w e can solve for 
^n+i recursively, starting with 



Assuming that ^> n is known, we determine ^ n +i/M from 
Eq. (|3.7[) which has a form similar to that of Eq. (|2.5[) . 
We use therefore the same method of Sec. [TT| to ob- 
tain V&n+i/M- This is repeated to obtain in succession 
^n+2/M, *n+3/M, ■ ■ ■ , *n+(M-i)/M, *n+i- Since the op- 
erators Kg commute, they can be applied in any order. 



IV. DISCUSSION OF ERRORS AND 
BOUNDARY CONDITIONS 



A. Errors 

In this section we discuss the errors as a function of the 
orders of the method, i.e., r and M. Let us separate the 
truncation errors due to the integration over space and 
those due to integration over time. At a given time t the 
spatial integration with the rth-order expansion yields a 
truncation error 



(3.7) 



2r 



(4.1) 
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where is assumed to be slowly varying with r. Ac- 
tually C*M = |V (2r) (arV)l/(2r!) for some x* in the range 
of spatial integration, and thus is model dependent. If 
we specify an acceptable error, the step size Aa; can be 
adjusted to obtain that error. Since Ax = (xq — Xj)/J, 
an adjustment of Aa; is equivalent to a change in J. Re- 
calling that xq — xj is fixed, we obtain 



Aa; 



x - xj 
J 



l/2r 



and 



constant 
J 2r ■ 



(4.2) 



(4.3) 



We have assumed that is approximately constant. 
The CPU time for the calculation is proportional to the 
number of basic computer operations in solving the ma- 
trix equation (|2.12|) . This involves elementary row oper- 
ations on r — 1 rows in J — 1 columns to bring the matrix 
to upper triangular form, plus J back substitutions to 
obtain the solution. Hence 



CPU time ex # operations ex Jr cx 



( e M)l/2r' 



(4.4) 



This form gives a minimum (optimum) CPU time which 
occurs when 



IneM 



(4.5) 



For the time integration we assume a truncation error 
independent of r. For a given r the error due to finite At 
has a first term in the expansion 



;(«3=CW(At) 



2A/+1 



(4.6) 



where again C^ M > is assumed to be a slowly varying func- 
tion of M. We note that the factor \ in the numerator 

and denominator of Eq. (]2.4[) is replaced by l/zi in 



each of the M factors 



of Eq. 



As M increases 



the average over different values of s of \zg '\, which wc 

denote as zi^g e , also increases. In fact zivgc is a linear 
function of M as is seen in Fig. [TJ The effective expan- 
sion parameter can be approximated by 2Ai/zavge rather 
than At and hence is proportional to At/M. Thus we can 
replace the relation of Eq. I]4.6[) by 



2JW+1 



(4.7) 



where the constant C^ M > is appropriately adjusted. If wc 



take the total time t. 



CPU time 



N At to be fixed, then 
1 

(M)\ 2M + 1 



(4.8) 



\z <M> \ 




FIG. 1: The average, the minimum and the maximum values 
of {|z<j M) |, s = 1 . . . M}, as a function of M. 



r, M increases from 1 through 5 or larger. For increas- 
ing M the CPU time continues to decline although the 
decrements become smaller at larger M. For increasing 
r there is a minimum depending on the specified error 
and beyond the minimum the curve shows a slow in- 
crease with increasing r. Superimposed on the curves 
are the CPU times (as dots) of a model calculation (see 
Sec. IV A]) , in which the numerical and exact solutions can 
be compared. Clearly the theoretical trends, including 
the minimum as a function of r, occur in the computed 
example. It should be noted that it "pays" to increase 
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FIG. 2: The normalized theoretical variation of the CPU time 
for a given error of 1.0 x 10~ 8 . The calculated CPU times for 
the example of Sec. IV Al are shown as dots. 



In Fig. [2] the curves of the (scaled) CPU times are plot- 
ted. Both curves clearly show the sharp decline when 



M indefinitely, whereas there is an optimum value of r 
which depends on the magnitude of eJ r \ 
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B. Boundary conditions 



EXAMPLES 



Below Eq. (|2.9j) we indicate that we set ip 



j,n — when 

j < or j > J, or when j goes out of range. These are 
appropriate boundary conditions when the wave function 
and its first r + 1 derivatives are zero at the boundaries, 
since it was assumed in the derivation of the method 
that all these derivatives exist. If however that is not 
the case, for instance, at the boundary of an (in)finitc 
square well or barrier where the second-order derivative 
does not exist, one must devise ways of incorporating the 
proper boundary conditions. 

One case of importance, which we discuss in the third 
example (see Sec. IV C| of the paper, is the case of radial 
behavior of a partial wave when angular momentum de- 
composition has been done. In the S'-wave case the wave 
function, defined only for nonnegative values of the radial 
coordinate, is zero at the origin but the first derivative 
is finite. Muller [22| discusses a related, but not identi- 
cal, situation. He considers three-point formulas for the 
Coulomb potential which lead to a radial wave function 
which is zero when the radial variable p = 0, but has first 
and second derivatives which are nonzero at p = 0. His 
approach can be adapted to the (2r + l)-point formula 
of this work. 

We treat this case by making the ansatz that the wave 
function behaves like an odd function about the origin 
and continues in the unphysical region of the negative 
radial variable. With this assumption we do not affect 
the behavior of the system at positive values of the radial 
variable, but all the required derivatives exist. Further- 
more, the wave function at negative j values can be com- 
bined with the ones with corresponding positive j values, 
so that the space need not be enlarged but instead the 
first few matrix elements of the matrix A can be changed 
to account for the boundary condition. This is achieved 
by replacing A by A' = A — B in Eq. (2.13) where 



B 



( a\ a 2 a 3 

a>2 0,3 Ct4 

CI3 0,4 CI5 

a r _i a r 

a r 





d r — 2 Q>r— 1 Or 

a r _i a r 
a r 













(4.9) 

A hard-core type potential could be dealt with in the 
same way. Different forms of boundary conditions are 
more complicated to implement, but Ref. [22j suggests 
an approach to including such boundary conditions. For 
the purpose of the radial wave function of a nonsingular 
potential, the above approach is sufficient. 



We consider three systems to which this numerical 
method may be applied. The first two allow us to make 
a comparison with the exact solution and to test the pre- 
cision of the numerical procedure. The third involves the 
time evolution of a quasi-stable quantum process. 



A. Oscillation of a coherent wavepacket 

The oscillation of a coherent state in the harmonic os- 
cillator well is described in Ref. [23| . The time evolution 
of such states has been discussed recently in connection 
with the quantum abacus. (See Ref. Q.) In that case 
there is a point interaction at the center of the oscillator 
well. It is of interest to consider narrow but finite-range 
interactions to simulate more realistic physical systems. 
In order to test the robustness of the quantum gates one 
needs a very stable numerical procedure. We test the 
precision of the numerical procedure by investigating the 
case without the central interaction, so that the numeri- 
cal results can be compared with the exact ones. 

The time-dependent Schrodingcr equation is 



. d . . 



2m dx 2 



l -Kx 2 ] 



(5.1) 



We consider the time evolution of the initially displaced 
ground-state wave function 



til 
rV4 



e -i" 2 (z-a) 2 



(5.2) 



where a 4 = mK/fi 2 , to = \J K/m, and a is the initial 
displacement. The closed expression for the time evolved 
wave function is 



til 
rV4 



exp 



~2^~ £ocosw<) 2 



— i(— + ££o sincjt — -£q sin2o;t) 



, (5.3) 



where £ = ax and £o = aa. We set h = m = 1, u) = 0.2, 
and a = 10. We choose our space such that x £ [xq, xj] = 
[-40, 40]. The period of oscillation is then T = 10tt. We 
allow the coherent state to oscillate for eleven periods 
before comparing the numerical solution to the exact one. 
The error is calculated as e-i using the formula [191 ] 



(e 2 ) 2 



dx \tp{x,ti) - '0o X act(^,tl)|' : 



(5.4) 



where t\ = 11T for our example. The results including 
the relative CPU time [34| are displayed in Table [TTT1 

In the above tests we have tried to obtain a precision 
better than 10 -8 . While varying the number of steps 
for the spatial integration, we kept the number of time 
factors per time step constant at 20. Given that the 
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TABLE III: Summary of computational time and errors in- 
curred by using the numerical integration procedure when the 
initial wave function is the displaced ground state. The last 
column indicates a relative CPU run time. The upper half of 
the table gives the effects of changing the number of spatial 
steps; the lower half the effects of changing the number of 
time steps. 



total space is fixed and spans 80 units, we adjusted the 
number of spatial steps J to give the required precision. 
We limited (arbitrarily) the maximum number of spatial 
steps to 2100. With M = 20 the 15-point formula (r = 7) 
is most efficient. When r < 4 (less than 9-point formula), 
we were unable to reach the precision criterion because 
of the imposed limit on J. It is clear from the trend 
however that the efficiency is significantly less for the 
lower r values. The 9-point formula is roughly half as 
efficient as the 15-point formula. 

The effect of different order time formulas as seen in 
the lower part of Tablc llIIl is even more dramatic. For the 
spatial integration we used the 21-point formula (r = 10), 
and varied the time-order formula, i.e., M, from 20 to 1. 
We see at least two orders of magnitude improvement in 
computational speed as M is increased over this range. 

A comparison with the standard CN approach (r = 
M = 1) is instructive. We considered the same sys- 
tem with xq = —25 and xj = 25, Ax = 0.005 and 
At = 0.5(Aa;) 2 . The standard CN method yielded an 
error of 62 = 7.1 x 10~ 5 when t = T/4 which increased 
exponentially to e 2 = 2.7 x 10~ 3 at t = 10T. Whereas 
the CPU time in Table IIIII is given in seconds, the CPU 
time required to complete this last calculation exceeded 
24 hours. 

The computed CPU times shown in Fig. [2] exceeded the 
"theoretic" values by increasing amounts as r increased 
beyond 10. This can be attributed to the approximate 
nature of the error analysis in which the model depen- 
dence of C*M (and C^) was neglected. In this exam- 
ple a more elaborate analysis could be done since the 
wave function is known analytically. In practical situa- 
tions where a numerical method is used the analytic wave 



function is usually not known and an estimate such as we 
have given here would be all that is available. The main 
point is that dramatic improvements result both theoret- 
ically and computationally when larger values of r and 
M are employed. 

B. Propagation of a wavepacket 

For this example we return to the work of Ref. § and 
consider the main features of that analysis with a view 
of determining the improvement brought about by the 
generalizations of this paper. This problem was revis- 
ited by Moyer [l8j] to illustrate the efficacy of the Nu- 
merov method and the use of transparent boundary con- 
ditions for the propagation of free-particle wavepackets. 
The authors of Ref. [8j consider wavepackets impinging 
on a square barrier and study their behavior in time. We 
consider first free wavepacket propagation (without po- 
tential), and second the reflection and transmission of a 
wavepacket by a smooth potential. 

Thus we first assume V(x) = and take as initial wave 
function 

1>{x,0) = (27T ( r2)- 1 /4 ei fco(a ; -xo) e -(.T-.To) 2 /(2ao) 2 i 

(5.5) 

(Note that our a is that of Ref. § divided by V2.) The 
wave function at later time is given by 

ip(x,t) = (27TCT 2 )- 1 / 4 [l + i^/(2m ( 7g)]- 1/2 x 

J — (x — Xo) 2 /(2er ) 2 + ik (x — x ) — ikk^t/ (2m) 1 
6XP 1 l + M/{2mal) J 

(5 r? 

We use parameters comparable to those of Ref. [8|. 
We set h = 1 and m = \. The coordinate range we take 
is from —0.5 to 1.5 rather than from to 1 since over 
the smaller space the normalization of the packet is not 
as precise as we require because the tails of the Gaus- 
sian are nonzero outside the (0,1) interval. We choose 
(To = 1/20, fco = 507r, At — 2(Ax) 2 , and allow as much 
time as it takes the packet to travel from Xq = 0.25 to 
around 0.75. For the final position the numerically cal- 
culated wave function is compared to the analytic one 
and e 2 is determined. In Table IIVI we list some of the 
computed results. We observe that the traditional CN 
method (Af = 1 and r = 1) has a low precision and that 
using greater J (smaller Ax) results in modest gain in 
precision. By using higher-order time formula one can 
make significant gain in precision (seven orders of mag- 
nitude) with no increase in computational time compared 
with that of Ref. Q. (Compare the first row to the last 
two rows of Table [IVJ) Rows 5 through 9 of Table |IV] 
illustrate the transition from less precise to more precise 
solutions. It is consistent with the finding of the authors 
of Ref. Q who use an M = 2, r — 3 method and obtain 
two orders of magnitude improvement of the results of 
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Refs. P, The results are sensitive to surprisingly high 
orders of Aa; and At. 

Another test using wavepacket scattering to show the 
efficacy of the higher-order approach is scattering from 
a potential. Rather than using the square barrier of 
Ref. [1] . w e consider the repulsive Poschl- Teller type po- 
tential [M HI of the form 



V(x) 



Ti 2 /3 2 A(A- 1) 
2m cosh 2 (3x 



(5.7) 



Since this potential does not have discontinuities the im- 
proved CN method works well with it. The transmission 
and reflection coefficients are known analytically. We 
can also compute them by considering the wavepacket 
Eq. (|5.5|) incident on the potential. Over a sufficiently 
long time the wavepacket will have interacted with the 
potential and transmitted and reflected packets emerge 
and travel away from the potential region. At that point 
we can calculate the probabilities of the particle repre- 
sented by the packet on the left and on the right of the 
potential; these probabilities correspond to the transmis- 
sion and reflection coefficients provided the packet is suf- 
ficiently narrow in momentum space. This means that 
one needs an initial packet which is wide in coordinate 
space. In our calculation we choose (3 = 1, A = 2.5, 
m = 1, and <7q = 10. This gives a spread in the inci- 
dent momentum-space wave packet of = 0.05. The 
width in momentum space of the reflected and transmit- 
ted wavepackcts also has this value. The domain of the 
x coordinates is from —300 to +300 and the initial posi- 
tion of the packet is at xq = —150 to ensure that there 
is no overlap of the initial packet and the potential. We 
find good agreement between the transmission and reflec- 
tions probabilities determined by plane-wave scattering 
approach and the time-dependent calculation as shown 
in Fig. G3 
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2000 
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2.12 


20 
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TABLE IV: Summary of computational parameters used to 
calculate the propagating free packet and compare it to the 
analytic wavepacket. 
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FIG. 3: The transmission and reflection coefficients as a func- 
tion of k when (3 = 1, A = 2.5, m = 1, and oo = 10. The sub- 
script "wp" indicates that the coefficients are obtained from 
the emerging wavepackets. The quantities without subscripts 
are calculated using the time-independent method. 



Long-time behavior of decay of quasi-stable 
system 



There are few analytically solvable models of the time 
evolution of unstable quantum systems 0, [2?]]. Real- 
istic systems need to be solved numerically. Long-time 
calculations are required for systems which require both 
nuclear and atomic time-scales such as ionization and 
bremsstrahlung due to radioactive decay of the nucleus 
of an atom [3, 0, i]. To study the short-time anoma- 
lous power law behavior of the decay and the long-time 
inverse-power law behavior the method of this paper is 
appropriate. This is especially relevant because of the 
recently observed violation of the exponential-decay law 
at long times [28| . 

To illustrate the numerical method discussed in this 
paper as applied to decaying systems let us consider a 
variant of the model with a i5-shell potential Q , but with 
the S function replaced by a gaussian. Thus in the S 
partial wave of a spherical system the potential is 



V{p) = 



A 



: exp \—{p — a) 2 jw 2 



(5.8) 



where p is the radial coordinate. This potential reduces 
to the J-shell potential, Vs(p) = XS(p — a) when w — > 
0. For small but finite values of w this potential leads 
to scattering results which are good approximations of 
those of the <5-shcll interaction [35| . Initially the quantum 
system is in the state 



ip(p,0) = ^/2/asin(7rp/a). 



(5.9) 



In our example we take 7i=l, A = 3, m=i, a = 1 and 
w = 0.10. Using the numerical method of this paper, 
including the modification of matrix A as described in 
Sec. HVBl to take care of the boundary conditions at p = 
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0, we determine the wave function at later times, i.e., 
tp(p,t). From that we obtain the nonescape probability, 
as a function of t, P(t) = \ip(p, t)\ 2 dp, which is shown 
in Fig. 2J It clearly shows the exponential decay-region 




2 4 6 8 10 12 14 



t 

FIG. 4: The nonescape probability as a function of time for 
the interaction with A = 3, a = 1, and w — 0.10. We also 
take h = 1 and m = i. 

in time, the inverse power-law behavior for long times, 
the deviation from exponential decay at short times, and 
the transition regions Q . 

The quadratic short-time behavior is seen in Fi g. [H 
as is the inverse power law behavior at long times [29j. 
Remarkably the decaying system can be studied in this 
manner for a time exceeding thirty half-lives. In Fig. O 
we plot the square of the absolute value of the wave 
function at times t = 5, 10, 15. Notice that the wave 




r 



FIG. 5: The square of the absolute value of the wave function 
as a function of p at times t = 5, 10, 15 for the same param- 
eters apply as in Fig. [4] The insert gives the t = graph as 
well as the scaled potential, Vi(p) — V(p)/8. 

function (packet) has three distinct regions: a precursor 
due to energy components of the initial wave function 
larger than that associated with the exponential decay, 



the main packet which corresponds to the exponential de- 
cay at the resonance energy, and the follower, which is a 
small blip that stretches in time (travels more slowly) and 
is due to energy components in the initial wave function 
which have lower energy that the resonance energy. If the 
maximum spatial coordinate, which is set at 800 for this 
calculation, were set at a smaller value, say 400, then 
one observes a fuzziness in the precursor of the right- 
most wave function. This can be attributed to the fi- 
nite space in which the wavepacket travels so that the 
fast precursor has been partially reflected from the right 
boundary and interferes with the wave front of the main 
wavepacket. The numerical parameters for this calcula- 
tion are Ap = 0.1, At = 0.02, r = 20 and M = 20. 



VI. REMARKS 

The generalized CN method that we have presented in 
this paper gives many orders of magnitude improvement 
in the precision of the results and several orders of mag- 
nitude in the computational time required to obtain the 
results. Clearly, since this method enhances the efficiency 
of the numerical calculations, it can be a significant tool 
for studying time-dependent processes. It goes beyond 
the improvement of Ref. Q in a systematic way. It also 
is an advance over the method of Ref. [l8[ , since the Nu- 
merov method has an error 0((Ax) 6 ), and it is difficult 
to see how it can be generalized systematically to higher 
order spatial errors. The generalized time evolution al- 
gorithm can be applied to Moyer's [H[ method. 

We have applied the method for r and M up to and in- 
cluding 20 for both. Having achieved significant improve- 
ments with these values of r and M we did not consider 
larger values although there does not seem to be a prac- 
tical reason that this cannot be done. Although the ap- 
proach seems to saturate at r around 10 (see Table HTT|) . 
there is no visible saturation in the time-evolution part of 
the problem. One would expect higher orders of spatial 
errors to be significant when the wave function and/or 
the potential has large spatial fluctuations. Even so, at 
r = 7 we are using a 15-point formula for the second- 
order spatial derivative, and it is surprising that a smaller 
number of points in the formula is not sufficient for op- 
timal results in this case. It should be noted that the 
higher order method as discussed in this paper are suit- 
able only for well-behaved, sufficiently diffcrcntiable so- 
lutions; these occur when the potential function is well 
behaved. As the authors of Ref. [ljj point out for sin- 
gular functions higher-order methods do not necessarily 
lead to greater accuracy. 

It should be noted that in this paper we consider pri- 
marily one-dimensional systems, but the method applies 
equally well to partial-wave equations of two- or three- 
dimensional systems. The study of the decaying quasi- 
stable state is an example of the latter. 

An interesting avenue to investigate further is the 
impact that this approach may have on two- or 
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three-dimensional systems, where the number of vari- 
ables involved is equal to the dimension. The 
Peaceman-Rachford-type approach [3(|, also known as 
the alternating-direction implicit method, of factoring 
the approximation of the time evolution operator may 
apply as it did in Refs. [3,[3l| or more recently in, for ex- 
ample, Refs. [H,[33j]. Using the Crank-Nicolson method 
the authors of Refs. [HI, [3l[ show that in two dimension 
the kinetic energy parts of the evolution operator factor- 
izes. Whether such factorization can be generalized in 
the spirit of the method of this paper is under investiga- 
tion. 

Calculations on one-dimensional multichannel systems 
indicate that this approach also leads to substantially 
greater efficiencies. 

Preliminar y st udies with the Numerov spatial integra- 
tion scheme [18[ and the generalized time evolution as 
described in this work indicate that significant improve- 
ments occur if one incorporates appropriate changes in 
the spatial step size for different regions of space. This 
is important in the case of discontinuous potentials and 
potentials that have great variation in some region and 
little or no variation in other regions. Furthermore it is 
well known that the wavepacket has large fluctuation in a 



(short-range) potential region and little variation in the 
asymptotic regions. A great savings in computational 
time can be achieved by using different space-step sizes 
in the different regions. One needs to investigate whether 
such variable step size can be incorporated in the gener- 
alized spatial integration scheme of this paper. The ap- 
proach that dealt with the discontinuous first- or second- 
order derivative in this paper and Ref. [22j is worth ex- 
ploring. We intend to study this in the future. 
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